Regulation of biofilm gene expression by DNA replication in Bacillus subtilis

Abstract Bacillus subtilis relies on biofilms for survival in harsh environments. Extracellular polymeric substance (EPS) is a crucial component of biofilms, yet the dynamics of EPS production in single cells remain elusive. To unveil the modulation of EPS synthesis, we built a minimal network model comprising the SinI‐SinR‐SlrR module, Spo0A, and EPS. Stochastic simulations revealed that antagonistic interplay between SinI and SinR enables EPS production in bursts. SlrR widens these bursts and increases their frequency by stabilizing SinR‐SlrR complexes and depleting free SinR. DNA replication and chromosomal positioning of key genes dictate pulsatile changes in the slrR:sinR gene dosage ratio (g r) and Spo0A‐P levels, each promoting EPS production in distinct phases of the cell cycle. As the cell cycle lengthens with nutrient stress, the duty cycle of g r pulsing decreases, whereas the amplitude of Spo0A‐P pulses elevates. This coordinated response facilitates keeping a constant proportion of EPS‐secreting cells within colonies across diverse nutrient conditions. Our results suggest that bacteria may ‘encode’ eps expression through strategic chromosomal organization. This work illuminates how stochastic protein interactions, gene copy number imbalance, and cell‐cycle dynamics orchestrate EPS synthesis, offering a deeper understanding of biofilm formation.

In this model, the copy number of kinA is gkinA = 2 because the kinA gene is proximal to the chromosome origin, while gi (i = kinA, spo0A, spo0B) equals 1 during DNA replication and 2 after DNA replication, i.e.,

Figure S4
. Process of burst identification.From top to bottom: raw data, initial burst data obtained using the find_peaks method, data after removal of erroneously identified small bursts, and data after merging bursts with the interval shorter than one cell cycle.
The reported statistical data are extracted from the bottom graph.
The statistical results on inter-burst interval (Tib) and burst duration (Db) in the main text were obtained using a consistent calculation approach.Initially, we applied a smoothing process to the time series.We then utilized the find_peaks function from the signal processing submodule of the built-in scientific computing library scipy in Python.
This method enables the identification of burst positions and burst widths.We further eliminated minor bursts falsely recognized by this method.Owing to protein dilution, we could not tell whether the number of molecules decreased between bursts of reporter fluorescence intensity occurring within a time window shorter than one cell cycle.Thus, we merged bursts occurring within a window shorter than one cell cycle.Ultimately, we acquired the resulting burst data.
When constructing histograms for noise analysis, we excluded all bursts with heights shorter than fluctuations in the background signal, i.e., the number of EPS proteins within a burst should exceed a threshold Nth (e.g., 10).Notably, the statistics of the burst duration, burst interval and proportion of EPS-expressing cells are insensitive to the concrete value of Nth provided it is small (see Figure S5).

Table S1. Parameters in the model
Parameters with µ = log2 h -1 and after division (i.e., the cell volume is equal to V0).Here Z represents EPS proteins.

Table S2. The update law for stochastic simulation
The simulation was performed using the Gillespie algorithm.
The number of proteins is halved per division.

Table S3. Genome length and chromosomal locations of sinI-sinR and slrR in 28 bacteria
The gene location is measured clockwise from the replication origin (0° = 360° = chromosome origin; 180° = chromosome terminus).Data were drawn from the PubMed Gene Database.Δθ with parentheses denotes the effective Δθ, which can be used to calculate tr.When two genes are positioned on the left and right sides of the chromosomal ring (Figure 1 C), their effective Δθ is not directly obtained by subtracting them.Instead, it is calculated by symmetrizing the two genes to the same side.The parameter values were adopted from Ref. (2).

Figure S1 .
Figure S1.Temporal evolution of gene copy numbers in individual cells.Time courses of gIR, grR and gr (= grR/gIR) (from top to bottom).The durations of the cell cycle and replication period are 3 h and 1.5 h, respectively.

Figure S2 .
Figure S2.Dynamics of the Spo0A-P level (S) and intermediate variable (m).The cell cycle period (Tcell) is set to 1, 3 or 6 h.The parameter values are also changed to make the amplitude of S(t) equal 0.8 or 1.2 times its default value (to differentiate between these cases, the Spo0A-P level is denoted by S′) .Of note, S and m are almost identical

Figure S3 .
Figure S3.Time courses of the numbers of network components on a single simulation trial.The label '(I+rR)-R' quantifies the relationship in size between the total amount of both SinI and SlrR and that of SinR, where a magnitude of 1 signifies a predominance of the former, and vice versa.Tcell = 1 h.

Figure S5 .
Figure S5.Burst statistics according to different baseline thresholds.The boundaries of EPS bursts are defined as the points where NEPS crossed a specified 'background' threshold (Nth).(A) The mean burst interval (Tib) and (B) burst duration (Db) and (C) the percentage (P) of EPS-expressing cells are influenced by the value of Nth.Because 1/Tib, 1/Db and P are low-pass filtered relative to Nth, the quantities above remain unchanged when Nth < 20.Here, Tcell = 1 h.

Figure S6 .
Figure S6.Variations in Spo0A-P dynamics under diverse conditions.Changing the angular difference (Δθ') between kinA and spo0F (A) or the cell cycle length (B) alters the amplitude of Spo0A-P oscillation.(C) The proportion (pm) of the time when m remains high (i.e., m > KIR) when the amplitude of S(t) is set to 0.8, 1, or 1.2 times its default value.The red solid line shows the duty cycle of gr pulsing versus Tcell when the amplitude of S(t) has its default value.

Figure S7 .
Figure S7.Heatmap for the average number of free SinR as functions of the complex formation rate constants.At Tcell = 1 h, the variations in the rate constant of SinI:SinR formation (kIR = 120, 300, 600 h -1 ) and that of SinR:SlrR formation (kRrR = 20, 50, 100 h -1 ) have minor impact on the time-averaged number of free SinR (NR).

Figure S8 .
Figure S8.Impact of the rate constants of complex formation on Tib, Db and P. Shown are the heatmaps for Tib, Db and P (from left to right) as functions of the rate constants of SinI:SinR formation (kIR) and SinR:SlrR formation (kRrR) for Tcell = 1 h (A), 3 h (B) or 6 h (C).kIR is greater than kRrR, ensuring that SinR has a higher binding affinity for SinI than SlrR.Overall, changes to kIR and kRrR have minor effects on Tib, Db and P.

Figure S9 .
Figure S9.Percentages of EPS-expressing cells in a colony under different conditions.

Figure S10 .
Figure S10.Impact of chromosomal arrangement of genes on EPS dynamics under diverse conditions.Mean values of Tib (A) and Db (B) and the percentage of EPSexpressing cells (C) versus the angular difference (Δθ) between slrR and sinI-sinR.Tcell = 1 h, 3 h, or 6 h, and the input S(t) takes the same form as in Figure 6B.
Z →  Z +  with probability